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ABSTRACT 

Cosmological tests based on cluster counts require accurate calibration of the space 
density of massive halos, but most calibrations to date have ignored complex gas 
physics associated with halo baryons. We explore the sensitivity of the halo mass 
function to baryon physics using two pairs of gas-dynamic simulations that are likely 
to bracket the true behavior. Each pair consists of a baseline model involving only 
gravity and shock heating, and a refined physics model aimed at reproducing the ob- 
served scaling of the hot, intracluster gas phase. One pair consists of billion-particle re- 
simulations of the original 500/i _1 Mpc Millennium Simulation of Springel et al.(2005), 
run with the SPH code Gadget-2 and using a refined physics treatment approximated 
by preheating (PH) at high redshift. The other pair are high-resolution simulations 
from the adaptive-mesh refinement code ART, for which the refined treatment includes 
cooling, star formation, and supernova feedback (CSF). We find that, although the 
mass funct ions of the gravity-o nly (GO) treatments are consistent with the recent cal- 
ibration of iTinker etatl (|2008l ). both pairs of simulations with refined baryon physics 
show significant deviations. Relative to the GO case, the masses of ~ 10 14 hr 1 Mq 
halos in the PH and CSF treatments are shifted by averages of —15 ± 1 percent and 
+ 12 ± 5 percent, respectively. These mass shifts cause ~ 30% deviations in number 
density relative to the Tinker function, significantly larger than the 5% statistical 
uncertainty of that calibration. 
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1 INTRODUCTION 



Deep cluster surveys offer the promise of tightly constraining 
cosmological parameters, including the nature of dark energy 
llHolder et ailhoOll; iLeyine et al1l2002l;jMaiumdar fc Mohrll2003t 
iLima fe Hull2004 12005| ; lYounger et altoofil; ISahlen et alJl200Sft . 



Realizing this promise requires accurate calibration of the ex- 
pected counts and clustering of massive halos, along with a careful 
treatment of how halo mass relates to the signals observed by such 
surveys. This logical division is reflected by two long-standing 
threads of effort, one focused on the emergence of massive struc- 
tures from gravity and the other focused on scaling relations of 
multiple signals within the population of massive halos. 

The fact that 17% of clustered matter in the uni- 
verse is baryonic ties these threads together. Non-gravitational 
physics i s required in massiv e halos, not simply to create 
galaxies (White & Rccs 1978) but also to reproduce scal- 
ing behavior of the hot, intracluster medium (ICM) ob- 
served in X-rays jEyrard fe Henrvl Il99ll; iBorgani etafl l200ll; 



Reiprich fc Bohringerl |2002| : IStanek et al.l l200rj : iNagai et al.l 
20071) . If a significant fraction of halo baryons become spatially 



segregated from the dark matter, either condensed within galax- 
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ies or disbursed from non-gravitational heating, then the gravita- 
tional development of massive structures will be altered, perhaps 
at the ~ 10% level, under strong baryonic effects. 

The spatial number density of halos, or mass function, ex- 
pected from Gaussian random initial conditions was originally de- 
rive d using a mix of analytic arguments and numerical simulations 
fe.g.lPress fe Schechterlll974l: iBond et ahlll99lt ISheth fc Tormenl 
Il999f) . Modern efforts focus o n providing fitting functions of in- 
creas i ng statistical preci sion i Jenkins et al.|| 200ll IWarren et al.l 
120061: ITinker et alj|200sl) . The recent ITinker et alJ l|200Sl) mass 
function (hereafter TMF), calibrated to a wide range of cosmo- 
logical simulati ons that include gas-dynamic, Maren ostrum simu- 
lations llYepes et alj|2007l ; [Cottlober fc YepesH2007f) , has pushed 
statistical errors to the level of 5%. 

To date, however, there have been few gas-dynamic sim- 
ulations that include a non-gravitational treatment of baryonic 
processes in volumes large enough to provide good statistics for 
hig h-mass halos. Calibr ation of the mass function at the level of 
the ITinker et al. using hydrodynamic simulations is too 

expensive to be feasible in the near-term. A less computationally 
expensive technique is to compare realizations o f fixed initial con - 
diti ons evolved with di fferent baryonic physics. I Jing et"al . (2006) 
and iRudd et al. (2008) employ this approach to study baryonic 
effects on the matter power spectrum, finding 2-10% modifica- 
tions of the matter power spectrum at scales k ~ lh Mpc - 1 . 
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Rudd etaj. | d2008h finds a halo mass function that is enhanced 
by ~ 10% relative to the dark-matter only case. Neither set of 
simulations were sufficiently large to properly probe rich cluster 
scales, however. 

In this letter, we take a similar approach to examine the 
effect of non-gravitational, baryonic physics on the cluster mass 
function. Specifically, we consider two pairs of gas-dynamic sim- 
ulations, each comprised of a treatment of the gas with gravity 
and shock heating only and a second, more com plicated treat- 
ment . One pair are Millennium Gas Simulations llHartlev et al.l 
l200Sl. MGS he reafterl . with SPH gas dynamics under Gadget-2 
(Springe! 2005), and the se cond pair are adaptive-mesh ART sim- 
ulations from Ruddl j2007l) . As the two simulations in each pair 
have the same initial conditions, we infer the effects of baryonic 
physics by comparing halos in the more detailed simulations with 
their gravity-only counterparts. The outline of this paper is as 
follows: in Section [2] we discuss the various simulations and their 
bulk cluster properties, comparing them to observations. Section 
[3] compares halo masses between corresponding halos in each pair 
of simulations, and discusses the effect on the total mass function. 
All halos masses are identified as M500, the mass of a spherical 
halo with radius rsoo and mean density 500p c (z), where p c (z) is 
the critical density of the universe. 



2 SIMULATIONS AND HALO SAMPLES 

We use two pairs of gas dynamic simulations, each pair run from 
a single set of initial conditions and differing only in the included 
physical processes. Our baseline simulations are evolved with only 
gravity and shock heating acting on the gas (hereafter GO) and 
a second has more complicated treatment of the gas physics dis- 
cussed below. 

Halos in all simulations are identified as spherical regions, 
centered on the peak of the dark matter distribution, where the 
mean enclosed density is 500p c (z). Our analysis focuses on 2 = 
and z = 1 samples, redshifts that roughly bracket the range im- 
portant for dark energy studies. When calculating the bulk clus- 
ter properties, we exclude from our sample halos that overlap 
with more massive neighbors; however, to remain consistent with 
iTinker et al halos that overlap but whose centers lie out- 

side their respective virial radii arc included in the mass function 
analysis. 



2.1 Millennium Gas Simulation 

The Millennium Gas Simulations (hereafte r MGS) are a pair o f 
resimulations of the Millennium Simulation iSpringel et al.l2005[) , 
a high-resolution, dark-matter-only simulation of a 500 h~ l Mpc 
cosmological volume. Like the original Millennium Simulation, 
the simulations were run with GADGET-2, which treat s the gas 
dyna mics with smoothed particle hydrodynamics (SPH) dSpringell 
2005). The MGS runs use a down-sampled version of the initial 
conditions of the Millennium simulation, with 5 X 10 s dark mat- 
ter particles, each of mass 1.422 X 10 10 h' 1 M Q , and 5 X 10 8 
SPH gas particles, each of mass 3.12 X 10 9 h _1 Mq. This mass 
resolution is about 20 times coarser than the original Millennium 
simulation, and the gravitational softening length of 25 h~ 1 kpc is 
correspondingly larger. The cosmological parameters match the 
Millennium Simulation: (fl m , Q\, h, n, erg) = (0.25, 0.045, 
0.75, 0.73, 1.0, 0.9). 

Complementing the aforementioned GO realization is an MGS 
simulation with cooling and preheating, denoted as PH . Pre- 
heating is a simple approximation that assumes high redshift 
galaxy formation feedback drove the proto-ICM gas to a fixed 
entropy level, after which the ICM evolve s under hierarchica l 
gravity tevrard fc Henrvlll99il : lKaiser|[l99ll : iBialek et al.ll20oil) . 



In our implementation, the entropy of each gas particle is in- 
stantaneously boosted to 200 keV cm 2 at z = 4. The gas is al- 
lo wed to radiatively co o l ther eafter using the cooling function 
of lSutherland fc Dopital Jl9Q3h . but the cold gas fraction is very 
small. The entropy level of the PH model is tuned to match bulk X- 
ray observations of clusters at redshift zero, as we discuss shortly. 

For both MGS models, we calculate bulk cluster properties 
with primary halos of mass M500 Ss 5 X 10 13 h- 1 M , yielding 
sample sizes of 2527 (PH) and 3446 (GO) at 2 = and of 475 (PH) 
and 818 (GO) at 2 = 1. 

2.2 ART Simulations 

Our second set of models are si mulated using the distributed - 
parallel hydrodynamic ART code l|Rudd et al.ll2008l : lRuddll2007l) . 
The simulations evolve a 240 3 h _3 Mpc 3 volume of a WMAP3- 
motivated cosmological model with parameters, (Q m , f^f,, S^Ai 
h, n, a 8 ) = (0.25, 0.042, 0.75, 0.73, 0.95, 0.8). The baseline GO 
simulation was performed using 512 3 dark matter particles with 
mass m v ~ 5.95 X 10 9 h~ 1 Mq and allow for 4 levels of refine- 
ment achieving a minimum cell size of 240 h' 1 Mpc/(512 X 2 4 ) Ri 
29 ft. -1 kpc. We then selected the 13 most massive halos in the 
simulation volume at 2 = 0, and resimulated at 1024 3 effec- 
tive resolution (m p fa 7.44 X 10 s h~ 1 Mq) the regions within 
5r v i r ~ 5 — 10 h~ 1 Mpc surrounding each cluster center. This 
simulation includes prescriptions for star formation a nd metal- 
dependent radiative cooling described in Rudd (2007]). For this 
simulation, the same 512 3 uniform mesh was used, but in the 
high-resolution regions 7 levels of refinement were used for a peak 
spatial resolution of ss 3.6 h~ 1 kpc. 

2.3 Baryon Census 

We begin by exploring the bulk properties of baryons in the mas- 
sive halo samples as a means of assessing the viability of the 
different physical treatments. 

Figure [T] shows baryon mass fractions within rsoo , normal- 
ized to the universal baryon fraction fij,/f2 m . In the GO simula- 
tions, all the baryons are in the hot intracluster medium phase, 
so that = /icm- I n t ne CSF simulation, gas is removed from 
the hot phase through radiative cooling and converted to stars. 
For these halos, we plot both the ICM mass fraction /icm an d 
the total baryon fraction, /(, = /icm + /condi where condensed 
baryons, / con di includes both stars and cold gas (T < 2 X 10 5 K). 
Although the PH model allows radiative cooling, the fraction of 
cold gas in our halo samples is very small, less than two percent 
of the baryons. 

The GO simulations display constant baryon fractions that 
are slightly suppressed from the universal mean value. The level 
of suppression is somewhat larger in the MSG halos compared 
to the ART sample, which is consistent with the difference be - 
tween SPH and grid codes reported by iKravtsov et al.l J2005h ■ 
The mass-limited samples have average baryon fractions at 2 = 
of f b = 0.89 ± 0.025 and 0.93 ± 0.038, respectively. The MGS 
baryon fractions are consistent with those in SPH simulations 
done at similar r esolu tion and measured within A c = 200 by 
ICrain et ah! 1120071 ) and lEttori et al.l ||2006| ). 

The baryon distributions in the PH and CSF simulations are 
more complicated. In the PH case, the entropy increase at 2 ~ 4 
causes the gas to expand, especially in lower-mass halos for which 
the characteristic entropy is lower, and raises the sound speed 
throughout the proto-ICM. The latter effect pushes the effective 
shock radius to large r valu es compared to the purely gravity- 
driven case jVoit et al.ll2003h . As a result, lower-mass halos retain 
a smaller fraction of their baryons within rsoo, leading to the 
mass-dependence seen in Figure [T] Since the characteristic halo 
entropy increases with time at fixed mass, the mean ICM gas 
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Figure 1. Baryon fractions as a function of total halo mass at 
z = 1 (left panel) and z = (right panel) are shown for the 
MGS halo samples in the GO (open circles) and PH (filled circles) 
treatments and for the ART samples in the GO (open triangles) 
and CSF (filled triangles) cases. For the latter, regular triangles 
show the ICM mass fraction while inverted triangles show the 
total baryon fraction (gas plus stars) within the halos. The shaded 
region shows the 90% confidence range of the mean, obs erved ICM 
mass fraction inferred for local, kT > 4keV clusters bv lVikhlininl 
(2006). For clarity, only a subset of the MGS samples are shown 
for masses below 4 X 10 14 /i -1 Mq. 



fraction at fixed mass increases from z = 1 to z = in the PH 
halos. At z = 0, the highest mass halos have baryon fractions 
suppressed by only 10% relative to the GO treatment. 

In the CSF halos, the hot gas fraction, /icmi i s comparable 
to that in the PH clusters of similar mass at Z = 1. However, the 
total baryon fraction in these halos is close to universal, due to the 
contribution of cold gas and stars. Unlike the PH models, the ICM 
mass fraction does not evolve with time, remaining approximately 
constant at ss 50% from z = 1 to z = 0, even as the clusters 
themselves grow by a factor of two in total mass. The total baryon 
fraction grows by 4% due primarily to the small increase in the 
ICM. The stellar component grows significantly from ~ 40% to 
~ 47% but is balanced by a corresponding decrease in the fraction 
of cold gas from ~ 10% to ~ 4%. 

As an empirical test of the models, we show in Figure [2] the 
scaling between bolometric luminosity Lbol an d spectral tempera- 
ture T sl for the 2 = halo samples of the MGS and ART-CSF sim- 
ulations. Th e models are c o mpar ed to a local sample of clusters 
compiled by lHartlev e t al. (2008|). As the local sample extends 
only to modest redshifts, z < 0.2, we do not apply evolutionary 
corrections to the observations. 

For the models, w e use the analytic approximation of 

iBartelmann fc Steinmet3 |l996) to compute L(, ; within r200 of 
each halo. For the GO and PH halos, we compute spectroscopi c 
temperatures, T s \, using the expression in lMazzotta et al] (2004). 
This expression is known to fail at low te mperatures , so, fo r the 
CSF clusters, we use instead the method of [Vikhlininl l |2006l ). Ad- 
ditionally, for the CSF clusters we exclude gas within 0.1r2oo an d 
within dark matter substruct ures to crude l y rep rodu ce the clump 
remov al procedure applied in lNagai et al] J2007h and|R asia et al] 
(2006). Applyin g this simple analys is proceedure to the simulated 
clusters used in lNagai et al. | (120071) give temperatures that differ 




1 10 
kT [keV] 



Figure 2. The Lbol — T^A relations for the MGS-G0 (open cir- 
les), MGS-PH (filled circles), and ART-CSF (filled triangles) sam- 
ples are compared to observations (small points) compiled by 
lHartlev et alj l|2008l l . The MGS halos are sub-sampled as in Fig- 
ure [T] 

by ~ 10% or less from the mock Chandra analysis. The measured 
X-ray quantities are sensitive to the choice of innermost radius, 
with larger cuts leading to simultaneously lower measured L DO i 
and T s i. 

Both of the non-gravitational physics models provide a bet- 
ter matc h to the observed da ta than the GO simulation. As dis- 
cussed in Hartley et al. (2008), the PH halos match the slope and 
normalization of the observed L — T relation well. The observed 
scatter is much larger, however, due primarily to the existence of 
cool cores in real clusters. 

The slope of the CSF halos also agrees with the observations, 
but the normalization and scatter are not well matched to the 
data. The normalization offset is partly due to the lower gas frac- 
tion seen in Figure p] but the spectral t emp eratures also play a 
role. As discussed inBorg ani et al ] l|2004h and lNaeai et al 1 l|2007l) . 
the temperature profile of the hot phase is steeper than observed 
in cluster cores, resulting in enhanced T s \ values. 

In summary, we have shown that both the PH and CSF simu- 
lations offer a reasonable match to the form of the L — T relation, 
but the overall baryon content of halos differs substantially be- 
tween the two. In the CSF simulation, star formation is overly 
efficient, so that nearly 50% of the baryons are in stars rather 
than in the hot phase. In the PH simulation, the stellar fraction 
is entirely neglected, but the net heating effect of early galaxy 
formation is assumed to be large enough to drive the halo baryon 
fraction substantially below the global value. Neither of these 
treatments is fully consistent with observations, but they repre- 
sent two extreme approximations for the true behavior. We next 
examine the effects that these treatments have on halo mass. 



3 HALO MASSES AND THE MASS 
FUNCTION 

Since both pairs of simulations are evolved from the same initial 
conditions, we are able to match halos between the realizations 
performed under the two physical treatments. In Figure [3] we 
show the fractional shift in mass that occurs under the PH and 
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Figure 3. Fractional mass difference in halo mass with respect 
to the GO realization. Circles show the mean shift for the PH halos 
at z = (filled) and z = 1 (open). Triangles show individual CSF 
halos at z = (filled) and 2 = 1 (open). The inset panel plots the 
cumulative radial mass difference for CSF halos (red) and PH halos 
(black) at 2 = (solid, with 1— a scatter) and 2 = 1 (dashed). 



CSF treatments, relative to the respective GO model, as a function 
of GO halo mass at redshifts 2 = and 1. The mean mass shift 
are plotted for MGS halos in mass bins. Individual clusters are 
plotted for the ART simulations at z = and 2 = 1. 

The PH halos experience a substantial decrease in mass rel- 
ative to the GO treatments. The magnitude of the 2 = frac- 
tional mass shift depends on halo mass, declining from 15% at 
10 14 h' 1 Mq to 5% at 10 15 h" 1 Mq. Although these mass shifts 
are mostly due to the change in gas fraction, there is also a dif- 
ference in dark matter structure that enhances the shift. All but 
the most massive ART-CSF halo show increased mass relative to 
the respective GO halos. At the mean mass of the sample at z = 0, 
3 X 10 14 h^ 1 Mq, the mean fractional mass shift is 0.117±0.015, 
including the outlier data point. Approximately 2% of this shift 
is due to the increase in baryon mass. The remainder is due to 
the change in halo structure br ought about by baryon cooling 
iGnedin et alj|20u4 lNagaill2006h . 

These mass shifts depend on the choice of scale, as shown 
in the inset of Figure \3\ For comparison with the ART— CSF ha- 
los we plot the mean mass profile for MGS halos in the range 
1 - 3 X 10 14 h- 1 Mq. Within the core, the mass difference be- 
tween matched halos in the MGS simulations is nearly 20%, but 
the mass difference approaches zero on scales significantly larger 
than T200- In the ART simulations, we also see that the mass 
shift is a strong function of scale: within the core it is very high, 
~ 80%, and approaches zero beyond T200- Because of this scale 
dependence, the magnitude of the mean mass shift and its evolu- 
tion with redshift is sensitive to our choice of A = 500p c (2). 

The shifts in mass seen with complex physical treatments 
will lead to changes in the mass function relative to the GO mod- 
els. For both MGS simulations and the ART-GO run, we compute 
binned space densities directly from the simulation counts. Fig- 
ure [4] plots these mass functions at redshifts 2 = and 2 = 1, 
and compares them to the TMF expectations for mean density 
contrasts equivalent to A c = 500, shown by the solid lines. To ac- 
count for differences in cosmology (primarily the difference in erg) 
we calculate the TMF for both cosmologies. From a fixed number 



Figure 4. The lower panel shows the halo mass functions for the 
PH (solid circles) and GO (open circles) versions of the MGS and the 
ART-GO model (open triangles) at redshifts 2 = (upper) and 
2 = 1 (lower). The solid black lines are the TMF expectations 
at these redshifts. Red bands shows the 90% confidence regions 
anticipated by the shifts in halo mass for the ART-CSF treatment. 
The panels above show the fractional difference in number counts 
between the measured mass functions and the TMF, at z = 
(top) and 2 = 1 (middle). We plot 90% jackknife errors for the 
ART-GO model, and Poisson errors for the MGS. Note that we 
have used the TMF for scaling the ART mass functions to match 
the MGS cosmology. 



density, we find the mass shift between the two cosmologies, and 
apply it to the ART-GO data for simple comparison with the MGS 
mass functions. 

The redshift zero GO mass functions match the TMF predic- 
tion quite well. The top panels show the fractional difference in 
counts between the simulations and the TMF, with the 90% sta- 
tistical calibration uncertainty of the latter shown by the solid, 
horizontal lines. We include 90% uncertainties on the data points: 
jackknife uncertainties as a measure of cosmic variance for the 
ART-GO sample, and Poisson uncertainties for the larger volume 
of the MGS simulation. The counts of both the ART and MGS 
models under GO treatment lie within the TMF expectations at 

2 = 0. 

At all masses, the PH halo mass function is suppressed with 
respect to the GO halo mass function, at a statistically significant 
level. At Af5oo ~ 10 14 h~ 1 Mq, the number density of PH halos is 
20% lower than the TMF prediction, a 4<x shift relative to the 5% 
TMF calibration error. At the very high mass end, ~ 10 15 M o , 
there is consistency with the TMF expectations. 

We do not have a complete mass function from the CSF simu- 
lation. However, we can anticipate the shift in halo number based 
on the mean shifts in halo mass presented above. Since the ART- 
GO models are consistent with the TMF, we derive CSF expec- 
tations by shifting the mass by fractional values given by the 
90% confidence range of the mean shifts shown in Fig [3] meaning 
AM/M = 0.126 ± 0.023 at z = 1 and 0.117 ± 0.024 at z = 0. 
We apply these shifts at mass scales probed by the CSF halos, 
M 50 o > 2 X 10 14 h -1 Mq. At 2 = 0, the positive shift in halo 
mass implies upward deviations in number density from the TMF 
expectation, at levels ranging from 10% to 60%. 
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4 DISCUSSION AND CONCLUSION 

Calibrations of the halo space density from ensembles of N-body 
and dissipationless gas dynamic simulatio ns now have very s mall 
statistical uncertainties, ~ 5% in number llTinker et al . 200^). At 
the high-mass end, this level of precision in number is equivalent 
to a precision in halo mass at the 2% level. Since baryons represent 
17% of the matter density, complex gas dynamics associated with 
galaxy formation physics could plausibly lead to effects on halo 
masses of more than a few percent. In this letter, we demonstrate 
that shifts approaching 10% in mass are possible, and that the 
sign of this effect is not yet understood. 

We use two extreme treatments of gas physics that are likely 
to bracket the range of behavior due to astrophysical processes in 
galaxy clusters. A simple assumption of preheating reduces the 
local baryon fraction in halos, thereby suppressing their mass at 
levels ranging from 15% at lO 14 ft- 1 M to 5% at 10 15 H^Mq. A 
more complete physics treatment with cooling and star formation 
increases the local baryon fraction and deepens the halo potential, 
thus enhancing halo mass, by an average of 12% at 10 14 ' 5 /i _1 M0. 
The effects of cooling and star formation on halo mass are qualita- 
tively consistent with the systematic enha ncement in small-scale 
powe r seen in previous simulations djine et al.ll2006t I Rudd et aU 
l2008t ). In both of the complex physical treatments we consider, 
the shifts in mass lead to statistically significant offsets in cluster 
counts from the TMF expectations. These shifts in mass depend 
on the choice of scale used in defining halos: in both treatments, 
the mass shifts are larger when identifying halos via higher den- 
sity contrasts. 

Although both the PH and CSF simulations provide fair 
matches to the mean observed L — T relation, implying the struc- 
ture of the hot gas phase is nearly correct, neither describes well 
the stellar content of clusters. The PH simulation ignores galaxies 
while the CSF simulation converts nearly ~ 50% of baryons into 
a large stellar component. Although it is tempting to dismiss the 
PH model due to its lack of detailed physics, a growing body of ob- 
servations, particularly t he ubiquity of stron g winds in moderate 
redshift DEEP2 galaxies llWeiner et al.ll2008T) and the remarkably 
simple evolution to z = 1.4 of the color of red sequence galaxie s 
seen in the Spiter /IRAC Shallow Survey l lEisenhardt et al.l2008t) , 
provide supporting evidence for a scenario in which the fireworks 
associated with galaxy formation in clusters is both rapid and 
effective. 

Improvements in the physical and computational modeling 
of cooling and star formation are needed to match the full set 
of observational constraints on the baryonic mass components of 
cluster halos. We have shown here that varying these treatments 
can affect total halo masses at levels up to ten percent. Improving 
the accuracy of the halo mass function calibration will therefore 
entail a suite of sophisticated gas dynamic simulations, not more 
or larger N-body simulations. 
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